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% METHOD OF FORMING A MODEL REPRESENTATIVE OF THE 
DISTRIBUTION OF A PHYSICAL QUANTITY IN AN UNDERGROUND ZONE, 
FREE OF THE EFFECT OF CORRELATED NOISES CONTAINED IN 

EXPLORATION DATA 

BACKGROUND OF THE INVENTION 

Field Of The Invention 

[0001] The present invention relates to a method of forming, from data 
obtained by exploration of a zone of a heterogeneous medium, a model 
representative of the distribution in the zone of a physical quantity (at least 
partly) free of the presence of correlated noises that may be contained in the 
data. 

[0002] The method applies A for example^ to the quantification of the acoustic 
impedance in an underground zone. 

BACKGROUND OF THE I NVENT I ON 

Description of the Prior Art 

[0003] The process cons i sting i n of seeking a model that adjusts to 
experimental measurements has been developed in nearly all the fields of the 
sciences or technology. Such an approach is known under various names: 
least-squares method for parameters estimation, for inverse problem solution. 
For a good presentation of this approach within the context of geosciences, one 
may for example refer to: 

Tarantola, A.: "Inverse Problem Theory: Method for Data Fitting and 
Model Parameter Estimation", Elsevier, Amsterdam, 1987. 
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[0004] It can be noted that the term "test- least squares" refers to the square 
of the norm in the data space for quantifying the difference between the 
response of a model (which is the image of the model by a previously selected 
modelling operator) and the data using , a cost function that has to be minimized 
to solve the problem. Using the square of the norm to define the cost function is 
just a practical convenience, and it is not fundamentally essential. Besides, 
many authors use, for various reasons, another definition of the cost function 
but this definition remains based on the use of the norm, or of a semi-norm, in 
the data space. Finally, w e hav e considerable latitude exists for selecting the 
norm (or the semi-norm) in the data space ( w e ar e in no way forc e d to the use of 
the Euclidean norm is not required ). In the case of noise-containing data, the 
solution can substantially depend on the choice made at this stage. For m More 
developments concerning this problem, one can for examp le r e f e r are referred to 
in the following publications: 

Tarantola, A.: "Inverse Problem Theory: Method for Data Fitting and 
Model Parameter Estimation", Elsevier, Amsterdam, 1987 ; Renard and Lailly, 
2001 ; Scales and Gersztenkorn, 1988 ; Al-Chalabu, 1992. 

[0005] The measuring results often contain errors. Modelling noises add 
further to these measuring errors when the experimentation is compared with 
modelling results: modellings are never perfect and therefore always 
correspond to a simplified view of reality. The noise will therefore be described 
hereafter as the superposition of: 

non correlated components (white noise for example), 
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correlated components, which means that the existence of a noise on a 
measuring sample translates into the existence of a noise of equal nature on 
certain neighbouring measuring points ; modelling noises typically belong to this 
category. 

[0006] When the data contain correlated noises, the quality of the model 
estimated by solving the inverse problem can be seriously affected thereby. As 
already mentioned, no modelling operator is perfect. It is therefore the work of 
the whole community of the people involved in the identification of parameters 
describing a model which is hindered by the existence of correlated noises. 
Among these people, seismic exploration practitioners are among the most 
concerned ones: in fact, their data have a poor or even very bad signal-to-noise 
ratio. This is the reason why correlated noise filtering techniques are an 
important part of seismic data processing softwares. The most conventional 
techniques use a transform (Fourier transform for example) where signal and 
noise are located in different areas of space, thus allowing separation of the 
signal and o ffrom the noise. For a A general presentation of the conventional 
methods for noise elimination on seismic data , on e may r e f e r to is discussed in 
the book by Yilmaz: 

Yilmaz, O. 1987: "Seismic data — Data proc e ss i ng Processinq ", 
Investigation in Geophysics No.2, Society of Exploration Geophysicists, Tulsa, 
1987. 

[0007] However, these techniques are not perfect: they presuppose that a 
transformation allowing complete separation of the signal and of the noise has 
been found. In this context, correlated noises are particularly bothersome 
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because they can be difficult to separate from the signal (which is also 
correlated) and can have high amplitudes. It is therefore often difficult to 
manage the following compromise: either the signal is preserved, but a large 
noise residue remains, or the noise is eliminated, but then the signal is 
distorted. These filtering techniques can be implemented before solving the 
inverse problem: they constitute then a preprocessing applied to the data. The 
quality of the inverse problem solution then widely depends on the ability of the 
filters to eliminate the noise without distorting the signal. 

[0008] The- This approach is introduced by Nemeth et al. in the following 
publication: 

Nemeth T. et al. (1999), "Least-Square Migration of Incomplete 
Reflection Data", Geophysics, 64, 208-221 which constitutes an important 
advance for the inversion of data disturbed by high-amplitude correlated noises. 
T hese The authors propose eliminating the noise by solving an inverse 
problem: after defining the space of the correlated noises as the image space of 
a vector space B (space of the noise-generating functions) by a linear operator 
T and the signal as the image space of a vector space M (models space) by a 
modelling operator F (assumed to be linear), they seek the signal in F(M) and 
the noise in T(B) whose sum is as close as possible (in the sense of the 
Euclidean norm or, on a continuous basis, of the norm L 2 ) to the measured 
datum. This technique constitutes an important advance insofar as it allows 
elimination (or at least very substantial reduction) of high-amplitude correlated 
noises (surface waves for example) that are difficult to separate from the signal 
by means of conventional filtering techniques. However, according to th e s e the 
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authors, an increase by an order of magnitude of the computing time required 
for solution of the conventional inverse problem, irer -that is without seeking the 
correlated component of the noise, is the price to pay for these performances. 
Furthermore, the result obtained with the method is extremely sensitive to any 
inaccuracy introduced upon definition of operator T: this is the ineluctable 
compensation for the high aptitude of the method to discriminate signal and 
noise. 

SUMMARY OF THE INVENTION 

[0009] The method according to the invention allows to estimate, from data 
obtained by exploration of an underground zone, a model representative of the 
distribution, in the zone, of at least one physical quantity, this model being 
a mp l y substantially f ree of the presence of correlated noises that may be 
contained in the data. The method essentially comprises the following stages: 

a) acquisition of measurements giving information about certain 
physical characteristics of the zone by following a predetermined experimental 
protocol, 

b) specification of a modelling operator which associates, with a 
model of each physical quantity, synthetic data that constitute the response of 
the model, the measurements and the synthetic data belonging to a data space, 

c) for each correlated noise referenced by a subscript j ranging from 
1 to J, selection of a modelling operator which associates a correlated noise 
with a noise-generating function belonging to a predetermined space of noise- 
generating functions, 
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d) specification of a norm or of a semi-norm in the data space, 

e) specification of a semi-norm in the space of the noise-generating 
functions spac e — for which each noise modelling operator establishes 
substantially an isometric relation between the space of the noise-generating 
functions space and the data space, 

f) definition of a cost function quantifying the difference between the 
measurements on the one hand and the superposition of the model response 
and of the correlated noises associated with the noise-generating functions on 
the other hand, and 

g) adjustment of the model and of the noise-generating functions by 
minimizing the cost function, by means of an algorithmic method taking 
advantage of the isometry properties of the noise modelling operators. 

[0010] According to a first implementation mode, the distribution as a 
function of the depth of the acoustic impedance in the medium is sought, the 
correlated noises affecting the data are tube waves identified each by 
parameters characterizing their propagation, the measured data are VSP type 
data obtained by means of pickups suited to detect the displacement of 
particles in the medium in response to a localized seismic excitation, the 
location of the pickups, the recording time and the time sampling points being 
defined, and the modelling operator selected associates the synthetic data with 
an acoustic impedance distribution as a function of the evaluated depth in 
traveltime and with the vertical stress measured as a function of time at the 
depth of the first pickup. 
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[0011] The cost function selected to quantify the difference between the 
measurements on the one hand and the superposition of the model response 
and of the correlated noises associated with the noise-generating functions on 
the other hand is for example the square of the semi-norm of this difference in 
the data space. 

[0012] According to an embodiment, adjustment of the model and of the 
noise-generating functions is for example obtained by means of a relaxation 
method for eliminating the unknowns corresponding to each correlated noise 
generating function, this relaxation method being implemented within the 
iterations of a quasi-Newtonian algorithm for calculation of the model. 

[0013] According to an embodiment, numerical calculation of the image of a 
model by this operator is carried out for example by numerical solution of the 1D 
waves equation for the model considered, by selecting values taken by the 
displacement of the particles at the locations of pickups and at the previously 
defined time sampling points, and by applying an operator likely to compensate 
for the spherical divergence and attenuation effects. 

[0014] According to an embodiment, the numerical noise modelling operator 
is a finite-difference centered numerical scheme for discretizing the noise 
transport equation, and the noise-generating function involved as the initial 
condition along the edge of said- the zone associates with each correlated noise 
a space of noise-generating functions space cons i sting of th o comprisinq 
support time functions in a given time interval. 
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[0015] Various examples of norms or semi-norms selected for the data 
space and the space of the noise-generating functions space are given in the 
detailed description hereunder. 

[0016] According to a second implementation mode, the distribution of 
disturbances, in relation to a previously selected reference medium, of the 
impedance and of the velocity in said- the zone of the medium is sought, the 
correlated noises affecting the data are due to multiple reflections whose 
kinematics and amplitude variations with the offset have been previously 
estimated, the measured data are picked up by seismic surface pickups, the 
location of said- the pickups, the seismic excitation mode, the recording time and 
the time sampling points being defined, and the modelling operator being 
defined via a linearization of the waves equation around the reference medium. 

[0017] The cost function selected to quantify the difference between the 
measurements on the one hand and the superposition of the model response 
and of the correlated noises associated with the noise-generating functions on 
the other hand is for example the square of the semi-norm of this difference in 
the data space. 

[0018] According to an embodiment, adjustment of the model and of the 
noise-generating functions is obtained by means of a block relaxation method 
for eliminating the unknowns corresponding to each correlated noise generating 
function, this relaxation method being implemented within the iterations of a 
conjugate gradient algorithm for calculation of the model. 
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BRIEF DESCRIPTION OF THE FIGURES DRAWINGS 

[0019] Other features and advantages of the method according to the 
invention will be clear from reading the description hereafter of an embodiment 
given by way of non limitative example, with reference to the accompanying 
drawings wherein: 

F i gur e Fig. 1 shows an example of data obtained by means of a vertical 
seismic prospecting method (VSP), 

F i gur e Fig. 2 shows a model of distribution of the acoustic impedance 
with depth, 

Figure Fig. 3 shows synthetic VSP data obtained on the basis of the 
impedance model of Fig.2, 

F i gur e Fig. 4 shows an example of VSP data contaminated by a single 
correlated noise, 

F i gur e Fig. 5 shows an example of VSP data contaminated by two 
correlated noises, 

F i gur e Fig. 6 shows the impedance distribution as a function of depth, 
obtained by inversion of the noise-containing data of Fig.4, 

F i gur e Fig. 7 shows the inversion residues (differences between the data 
of Fig.4 and the seismic response of the model of Fig.6), 

F i gur e Fig. 8 shows the impedance distribution as a function of depth, 
obtained by inversion of the noise-containing data of Fig. 5, 
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Figur e Fig. 9 shows the corresponding inversion residues (differences 
between the data of Fig. 5 and the seismic response of the model of Fig.8), 

Figuro Fig. 10 shows the impedance distribution as a function of depth, 
obtained by inversion of the noise-containing data of Fig.4 by seeking the 
correlated noise in a_form of the superposition of two correlated noises having 
inaccurate propagation properties, frer that is different from those of the noise 
appearing in the data of Fig.4, 

Figur o Fig. 11 shows the corresponding inversion residues (differences 
between the data of Fig.4 and the seismic response of the model of Fig. 10), 

Figur e s Figs. 12A~ and 12B respectively show seismic data with multiple 
reflections and the seismic response of the model obtained after conventional 
linearized inversion, 

Figur e s Figs. 13Ar- and 13B respectively show an example of impedance 
distribution and of propagation velocity for which the seismic data of Fig.12A 
constitute the seismic response, and 

Figures Figs. 14A— and 14B respectively show the comparison of the 
seismic responses of the models obtained by conventional inversion ( Fig. 14A) 
and by inversion according to the proposed method involving picking of the 
moveout of the multiple and an estimation of the amplitude variations with the 
offset (Fig. 14B). 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
[0020] Formulation of the problem 
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[0021] W e hav e d Data is available resulting from a sampled measurement of 
a function. This function depends on several variables (space or space-time 
variables for example). 

[0022] These data, called d, correspond to measurements carried out to 
obtain information about a model m. They contain various noises: 

correlated noises (or a superposition of correlated noises), 

non correlated noises. 
[0023] The problem is to determine quantitatively model m (or functions of 
this model) from datum data d. 

[0024] We therefore select a A modelling operator F is selected (linear or not) 
which associates with model m the response of the model. This operator 
actually models the real physical phenomenon only imperfectly. This is 
essentially (but not exclusively) the reason why correlated noises appear in the 
data: these noises correspond to a signal related to the model but the relation 
between them appears to be too complex to be included in modelling operator 
F, which we want to should keep relatively simple for various reasons. 

[0025] Under such conditions, w e mod elthe noise is modeled by introducing 
one or more noise-generating spaces Bj (j ranging from 1 to J) and the 
associated noise modelling operators called T,. In order to find model m from 
data d, w e propos e seeking this model is proposed as the solution to the 
optimization problem. 

mill j C(m 9 fi l ,p 2 ,...fi J ) = ^ (1) 

{m,fi l ,fi 2 ,...fij)eMxllB J » M 
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where || \\ D is a norm (or possibly a semi-norm) in the data space. 

[0026] The above formalism relationship differs from the approach proposed 
by Nemeth et al. (1999) in that: 

modelling operator F is not necessarily linear, 

w e ar e justifi e d in considering a superposition of correlated noises of 
different types is justified (in the sense that they are modelled by different 
operators referenced by subscripts j), at the price of a marked increase in the 
number of calculations concerning the number of unknowns (which particularly 
complicates the solution of the optimization problem) as well as the number of 
noise modellings to be carried out, 

w e r e serv e the possibility of selecting norm || \\ D or possibly the semi- 
norm to suit our convenience in the data space is reserved . 

[0027] The choice that w e a ro go i ng to mako is made for || \\ D will be 

imposed by considerations linked with the efficiency of the solution method 
(described below), this efficiency allowing to— dealing with the case of a 
correlated noise superposition. This is an important opening: besides the 
possibility of processing data containing correlated noises with complex 
characteristics, w e can in this way ov e rcom e t he difficulties are overcome which 
are encountered in Nemeth et al.'s approach, linked with the high sensitivity of 
the result to an inaccuracy introduced upon definition of the noise modelling 
operator— To te overcome this difficulty, w e just h a v e to s eo k the noise is 
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sought by superposition of noises associated with approximate modelling 
operators Tj. 

[0028] It can be noted that selection of the square of norm || \\ D involved in 

the optimization problem is not essential: we do not chang e t he solution is not 
changed by selecting instead a function which, like the square function, is an 
increasing function on 5R + (or a decreasing function provided that the min is 
changed into the max). 

[0029] The solution method 

[0030] The method consists i n suitably s ele ct i ng selects : 
norm (or semi-norm) || \\ D in the data space, 

norms || || in each space Bj 

so as to obtain a high-performance algorithmic solution to optimization problem 
(1). The expression "suitably select" means see to it that each operator T, 
constitutes an isometry (or an isometry approximation) for norms || || and | \\ D 

respectively in the initial space and in the end space. The main idea in the 
determination of norms a allowing to mak e allows making each operator Tj to be 
isometric, is based on the existence of energy conservation type equalities 
verified by the solutions to the partial differential equations (or of discrete 
energy equalities for the solutions to certain numerical schemes used for 
discretization of these equations). 
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[0031] If w e can m a k e such choices are made , it is possible to make the cost 
function minimization very simple by eliminating the unknowns associated with 
the correlated noise generating functions (determination of Schur's 
complement). This elimination is peformed by means of a suitable algorithm 
such as the block relaxation method, a block being associated with a noise- 
generating function, or the conjugate gradient method with symmetrical-block 
Gauss-Seidel preconditioning, etc., all of them algorithmic implementation 
methods well-known to tho man sk i l le d in the art, presented for example by: 

Golub, G.H. et al. (1983) "Resolution Numerique des Grands Systems 
Lineaires" (Eyrolles). 

[0032] Simplified minimization of the cost function allows to consid e rably 
fedwe considerable reduction of the time required for numerical solution. 

[0033] First implementation example applied to 1D inversion of VSP data 
contaminated by a tube wave 

[0034] VSP (Vertical Seismic Profile) data belong to conventional well survey 
acquisitions. Tbe^ VSP data have multiple applications, one of the main ones 
being their use for establishing a connection between seismic surface data and 
logging measurements. The inversion of VSP data as described for example by 
Mace, D. and Lailly, P. (1984) A Solution of an Inverse probl e m Problem with 
1D Wave Equation appli e d Applied t o the Inversion of Vertical Seismic Profile ; 
Proc. of the 6 th Intern. Conf. In "Analysis and Optimisation of Systems", Nice, 2, 
309-323 has been developed to that end. In the 1D inverse problem, we 
assume that the subsurface model is assumed to be laterally invariant and that 
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a plane-wave excitation propagates vertically. The unknowns of the problem are 
the acoustic impedance distribution as a function of depth (measured in vertical 
traveltime and over an interval starting at the depth of the first receiver) and the 
seismic excitation mode (precisely the time function characterizing the boundary 
condition at the depth of the first pickup). This is a non linear inverse problem, 
the VSP data depending non linearly on the impedance distribution. 

[0035] VSP data are often contaminated by a correlated noise of very high 
amplitude4-4fre from a fluid wave. This wave is guided in the mud that invades 
the well. Its propagation velocity is low in relation to the propagation velocity in 
the rocks surrounding the well. Furthermore, the propagation velocity of this-the 
fluid wave depends on the frequency (dispersive propagation). Finally, tWs-the 
fluid wave reflects notably at the well bottom and +t-furthermore gives rise to 
another propagation mode. Typical VSP data are given in the figure below: the 
main characteristics of the tube wave can be seen. Its amplitude and the 
extension of the contaminated zone make it difficult to use the signal contained 
in the many seismic samples contaminated by the fluid wave. 

[0036] W o d e scribe h Hereafter the implementation of the method for 1D 
inversion of the-VSP data is described . 

[0037] Experimental protocol 

[0038] It consists i n spoc i fyinq specifies : 

the various depths at which the receivers are positioned in the well 
(assumed to be vertical): in otH^experiments, the pickups cover the range of 
depths measured in one-way time [x m j n = 0.05s ; Xmax = 0.2s], with a receiver 
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every Ax=1 ms (in one-way time) ; we-thus ofetaifhsamples numbered from 0 to 
1 are obtained , 

the physical nature of the measurement performed by these- the pickups: 
in ew-experiments, the pickups measure the vertical displacement resulting 
from the wave propagation, 

the recording time: in ou^experiments, the pickups measure the 
vibrational state over the time interval [0, T=1.5s], these data being sampled 
every At=1 ms. We4Thus obtain samples numbered from 0 to N are obtained . 

[0039] We c a ll The data dj n jsjhe data sample recorded at the depth Xmm + i 
Ax and at the time n At. We have o Of course: T = N At and x max = x min + I Ax. 

[0040] Acquisition of data in accordance with tbis- the experimental protocol 

[0041] Our o Experiments were carried out from synthetic data, ^ that is 
obtained by numerical modelling, itself carried out in two stages: 

[0042] Creation of synthetic data containing no noise: this creation was 
performed by numerical solution of the 1D acoustic wave equation, the model 
being the acoustic impedance distribution shown in Fig.2 (the depth is 
measured in one-way time) and the excitation mode (Neumann's boundary 
condition at the depth 0) being a conventional Ricker wavelet. 

[0043] The synthetic VSP data resulting from this modelling are shown in 
F i gur e Fig. 3: the vertical axis representing the observation time is graduated 
from 0 to 1.5 s and the horizontal axis representing the depth (in one-way 
traveltime) of the various pickups is graduated from 0.05 to 0.2 s. 
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[0044] Addition of noise to the data 

[0045] We a re go i ng to disturb t The VSP data are disturbed thus calculated 
by adding thereto, on the one hand, a random noise (but filtered in the seismic 
frequency band) of rather high amplitude and, on the other hand, one or more 
correlated noises of very high amplitude. Each noise will propagate downwards 
with a low velocity which depends on the frequency (dispersive propagation): 
these correlated noises are supposed to represent fluid waves. They have been 
modelled by numerical solution (use of the finite-difference centered scheme) of 
a transport equation with a constant propagation velocity. Dispersion of the 
waves has been obtained using large discretization intervals in the finite- 
difference scheme. In the case of the superposition of several correlated noises, 
the propagation velocities associated with each noise are different. Figures 4 T 
and 5 show the VSP data contaminated by a single correlated noise (Fig.4) and 
by the superposition of two correlated noises (Fig. 5). 

[0046] Specification of the modelling operator 

[0047] Modelling is carried out by solving the wave equation: 

with Neumann's boundary condition: 

-oix^^-{x^i)=Kt)on[0j} 
ax 

and the initial conditions: 
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y(x,0) = (x 9 0) = 0on [x^ ,+qo] 
dt 

[0048] Definition of the modelling operator first requires definition of the 
observation operator which formalizes the experimental protocol described in 
§1 . The observation operator thus is the operator which associates with function 
y(x,t), solution to the above system, samples y(Xj,t n ) for Xi=x min + i Ax and t n =n 
At, i=0,...,l ; n=0,...,N. 

[0049] The modelling operator th a t w e cal l identified as F(m) thus is the (non 
linear) operator which associates with function pair m=(c(x),h(t)) samples y(Xi,t n ) 
for Xj=Xmin + i Ax and t n =n At, i=0,...,l ; n=0,...,N. 

[0050] Specification of the modelling operator associated with each 
correlated noise 

[0051] We specify in th i s paragraph t The procedure is specified herein that 
can be used for modelling each correlated noise. Each correlated noise will be 
characterized by correlation directions specific thereto, which w e assume are 
assumed to be known, at least approximately: these correlation directions are 
specified by means of a field of correlation vectors Cj(x 9 t) of components 

(Cj x (x,t), Cj^x.t)) which has to be defined at any point of the domain [x min , 
Xmax]X[0,T]. W e can, iln an equivalent way, sp e c i fy the correlation directions can 
be specified by means of a penc i l of ^correlation lines that will represent the field 
lines associated with c y (*,*)■ Such correlation lines can be obtained according 

to the technique described in pat e nts EP- 35 4 .1 12 ( US- Patent 4,972,383^- filed 
by the applicant, from a pick of several phases characteristic of the noise, the 
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information being extended to the whole domain [Xmin, XmaxlXfO,"!*] by means of 
conventional interpolation and/or extrapolation procedures. Specification of the 
correlation directions amounts to specifying the propagation velocity distribution 
of the noise Cj(x,t) that is related to the correlation vector by the relation 


wave amplitude variations during propagation of the wave (a different situation 
is presented for the 2 nd application example). Modelling of a correlated noise as 
specified above is based on the following observation: a wave propagating 
along correlation lines without amplitude change is solution to the transport 
equation: 


[0052] In this context, we-introduce d are : 

a space Bj of noise-generating functions which will be the space of 
support functions pj(t) in [tj min ,tj max ] <= [0,T] and that w o o xtend is extended by 0 in 
the whole interval [0,T], 

noise modelling operator Tj: 

TjiPjiOeBj-^bjiXtOeD 

solution to the transport equation 

Vbj (x 9 t),cj (x, t) = 0 in [x^ 9XmaK ] x [0, T] 

and meeting the initial conditions: 


c j (x,t) = 



c;(x,t>- 



Here it is assumed that there are no 


^h(x, t) . c{x, t) = 0 
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IZ>.(jc,f = 0) = 0 

[b (x = x^t) = Pj{t) 

[0053] It can be noted that the geometry of the pencil of correlation lines (or, 
which amounts to the same thing, the propagation velocity distribution) cannot 
be any geometry: the correlation lines cannot intersect, or the transport 
equation solution problem would not be correctly set. For the same reasons, the 
correlation lines cannot intercept twice the edges on which the initial conditions 
are specified. The latter consideration can lead (in order to model a correlated 
noise corresponding to an upgoing wave) to specify the initial condition no 
longer at x=x min but at x=x max or even at an intermediate value, even if it means 
decomposing Cauchy's problem into two problems associated with the domains 
located on either side of this intermediate edge. 

[0054] We- lt is furthermore assumed that the geometry of the correlation 
lines is such that (more general situations can be considered, but the following 
results then take a different form): 

cj(a?,*) = fc(ar) x r(t) 
ou encore : t) = ^(x) est c^x.t) = ^ 

[0055] Finally, we -it is assumed that the pencil of correlation lines that 
intercept at x=x min interval [tj min ,tj max ] do not intercept at t=T interval [x min ,x ma x]. 

[0056] In fact, the transport equation is solved by means of a numerical 
scheme. It is for example possible to use the conventional finite-difference 
centered scheme. We ^ Therefore soloct discretization intervals Ax/ and At/ are 
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selected (which are not necessarily equal to Ax and At, but which we-are 
assumed to be submultiples of these quantities, even if more general situations 
can be considered) and w e i ntroduce a grid is introduced whose nodes are the 
points of coordinates (x min + i'Axj , n'At/) with i' and n' e N. The conventional 
finite-difference centered scheme explained hereunder allows, from the initial 
conditions, to calculate step by step the different values taken by function bj(x,t) 
at the various nodes of the grid (to simplify the notations, w e h a v e left out 
subscript j associated with the correlated noise considered is omitted ). 

2^T [*V+1 + b i> ~ b i'+l ~ b i' J + 

[0057] In the above formula, quantity ar n represents the evaluation of 
quantity a (a is a generic term) at the point of coordinates (x mjn + TAXj , n'At/). 

[0058] Of course, the use of discretization intervals Ax/ and At/ sufficiently 
small in relation to the wavelength is required if w e want the numerical scheme 
is to give a precise approximation of the function solution to the transport 
equation. 

[0059] However, the use of larger discretization intervals (but still smaller 
than the wavelength) allows to — modelling more complex propagation 
phenomena by making the wave propagation velocity dependent on the 
frequency: dispersive propagations such as the propagation of tube waves can 
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thus be modelled. Selection of suitable discretization intervals for modelling a 
particular mode of the tube wave can be made: 

[0060] By determining, after examining the data, the propagation velocity of 
this mode (which depends here not only on space and time, but also on the 
frequency) ; tomography techniques can be used to that end, such as those 
presented in the following publication: 

Ernst, F. et al. (2000) ; Tomography of Dispersive Media ; J. Acoustic 
Soc. Am., 108, 105-116. 

[0061] By knowing the dispersive propagation properties associated with the 
numerical scheme, properties that follow from an analysis of the numerical 
dispersion relation, as described for example in the following publication: 

Alford, R.M. et al. (1974), Accuracy of Finite Difference Modeling of the 
Acoustic Wave Equation ; Geophysics, 6, 834-842. 

[0062] Specification of a norm or of a semi-norm in the data space 

[0063] A first choice consists in tak i nq takes as the norm in the data space 
any discrete norm which is meant to be an approximation (via a quadrature 
formula) to the norm defined on a continuous basis by: 
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[0064] We can also us e tT he semi-norm, which is a better choice as w e shall 
seediscussed in the next paragraph may also be used : 



where u is any vector of the data space (i and n are the indices representing the 
sample numbers respectively in space and in time). 

[0065] Besides, other choices are also possible. 

[0066] Specification, in each space Bj of noise-generating functions, of a 
norm for which each noise modelling operator 7} establishes an isometric 
relation (or the approximation of an isometric relation) between the noise- 
generating functions space and the data space (provided with the semi-norm as 
defined in §5) 

[0067] W e sp e c i fy i n this par a graph t The procedure herein specifies defining 
th a t can be used to defin e the norm in the noise-generating functions space 
associated with each correlated noise. This procedure being the same for each 
noise, w o d o scr i be i t gener i ca ll v the description is generic and leaves out, in the 
formulas, subscript j associated with the noise considered. 

[0068] A first choice consists i n prov i ding provides space of each noise- 
generating functions space B with any discrete norm meant to be an 
approximation (via a quadrature formula) to the norm defined on a continuous 
basis by: 
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in which case, using the fact that, by means of the hypotheses set out in §4, we 
hav e bj(x,T)=0 for x e [x m i n ,x ma x], linear operator Tj performs an isometry 
between Bj and D. However, since calculation of TjPj is carried out numerically, 
the isometric character of operator T will not be perfect but only approximate, 
the approximation being all the better as discretization intervals Ax, At, Ax/, At/ 
are small. The fact that w e then hav e only approximately an isometry wiH-is 
available leads to lower performances concerning the algorithmic 
implementation presented in p a r a graph 8 hereafter. 

[0069] A better choice cons i sts in providinq provides each noise-generating 
functions functional space Bj with the semi-norm: 



[0070] In this case, operator Tj rigorously performs an isometry between Bj 
and D, at least if u { N = OVi = 0,1-1 . This results from the discrete energy equality: 

'E-4rk' +»o) 2 ]=EZ-^[("f +1 +<) 2 ] 

n=0 * /=0 n=0 ' 

^ At K M *" +M * J J A t 


/=0 *=0 
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[0071] an equality valid in the case (simplified for presentation's sake) where 
Ax* = Ax and At' = At and if £.1/2 = £1/2 has been defined. 

[0072] Definition of the cost function 

[0073] We dofin e t The cost function is defined by the formula as follows: 

C(m, fi x ,fi 29 .. ,/3j ) = + £ Tjfij - <f j 

[0074] Seeking the model and the noise-generating functions minimizing the 
cost function, this search being carried out by an algorithmic method taking 
advantage, explicitly or implicitly, of the isometry properties of the noise 
modelling operators (see 6) 

[0075] The algorithmic method uses the specificity of Schur's complement, a 
specificity associated with the isometric character of operators Tj. This can be 
done for example by realizing that minimizing C(m,Pi, P2,..., pj) amounts to 
minimizing C(m) = C{m,p x {m\P 2 {m)^..Pj{m)) where (p l (m),P 2 (m),...p j (m)) 

minimizes C(m,Pi, P2,..., pj) for given m. It can be noted that the Hessian 
associated with this quadratic form is Schur's complement C'(m) can be 
minimized by means of any optimization method such as the BFGS version of a 
quasi-Newtonian method. 

[0076] Determination of (P(m) 9 p 2 (m),...Pj(m) is very easy thanks to as a 
result of the isometric character of operators Tj, notably if choices (2) and (3) 
have been respectively selected to define || and || || . If there is only one type 

of noise (J=1) and if the isometry is perfect, determination of p x (m) simply 
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requires evaluation of the gradient of the quadratic form. If there is a 
superposition of several noises (J>1), various algorithms can be considered to 
take advantage of the isometric character of operators Tj. (f? l (m) 9 /3 2 (m) 9 .../3 J (m) 
can for example be determined by means of a block relaxation method (Gauss- 
Seidel alias), a block being associated with a group of unknowns /?-(/w). a 

relaxation iteration simply requiring (still in the case of a perfect isometry) a 
single evaluation of the gradient of the quadratic form associated with the block 
considered. For problems made difficult by the proximity of the correlation 
directions of the various noise types, a preconditioned conjugate gradient 
method can be considered, the preconditioning matrix being of symmetrical 
block (relaxation alias) Gauss-Seidel type. Then again, preconditioning is very 
easily implemented— and solution of the problem associated with a block only 
requiring evaluation of the gradient of the quadratic form associated with the 
block considered. 

[0077] Figures Figs. 6 to 1 1 show the results obtained with three types of 
experiments using the method. 

[0078] In the case of an inversion of noise-containing data (Fig.4) both by 
superposition of a correlated noise and of a random noise, Fig. 6 gives the 
impedance distribution obtained, and Fig. 7 shows the inversion residues. 

[0079] In the case of an inversion of noise-containing data (Fig. 5) both by 
superposition of two correlated noises and of a random noise, Fig. 8 gives the 
impedance distribution obtained, and Fig. 9 shows the inversion residues. 
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[0080] In the case of an inversion of noise-containing data (Fig.4) both by 
superposition of a correlated noise and of a random noise where the 
propagation properties of the fluid wave are not precisely known, it is natural to 
seek a correlated noise model cons i sting of th e sup e rpos i t i on which superposes 
o^two correlated noises having inaccurate propagation velocities but by which 
the true propagation velocity of the fluid wave is framed. Fig. 10 gives the 
impedance distribution obtained T and Fig.1 1 shows the inversion residues. 

[0081] Second example of application of the method to determination of a 
model by linearized inversion of multi-offset data contaminated by multiple 
reflections 

[0082] This application differs from the first one mainly in that it illustrates the 
possibilities afforded by the method for taking into account noise amplitude 
variations along the correlation directions. Another difference lies in the 
propagation properties of the correlated noises, propagation of the multiple 
reflections being not dispersive. 

[0083] Linearized inversion is one of the sophisticated methods used by 
geophysicists to obtain a quantitative model of subsurface heterogeneities, this 
type of information being p re c ious — important notably for reservoir 
characterization. The data given here are surface seismic data. This method is 
described for example in the following publication: 

Bourgeois, A. et al. (1989) "The Linearized Seismic Inverse Problem: aft 
An Attractive Method for a Sharp Description of Strategic Traps": 59 th Ann. 
Intern. Mtg. Soc. Expl. Geoph. Expanded Abstract, pp.973-976. 
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[0084] Furthermore, it is necessary to have a reference model, consisting (in 
the context of an acoustic model) of a velocity distribution and of an acoustic 
impedance distribution: it is the model around which a linearization will be 
carried out in the waves equation (Born's approximation). In order to simplify the 
presentation, w e take the 1D case is taken but this hypothesis is all the less 
essential as the practical applications mainly concern 3D seismic methods. In 
the 1D context, all of the seismic information is contained in a single shotpoint. 
Besides, the reference model is also one-dimensional: it only depends on the 
depth. The unknowns of the problem are the acoustic impedance disturbance 
and velocity distribution as a function of depth. On the other hand, unlike the 
first application example, the seismic excitation mode (and notably the seismic 
wavelet) are assumed to be known. It is an inverse problem made linear by 
means of linearization. 

[0085] Surface seismic data are often contaminated by multiple reflections 
which may have a rather high amplitude. Their kinematics, which is 
characterized by moveout in the 1D context, is generally quite different from that 
of primary reflexions: these two types of waves have generally travelled through 
geologic layers having different propagation velocities. This is the reason why 
linearized inversion is in itself an efficient technique for attenuating multiple 
reflections: the kinematics of primary reflections (the events modelled by Born's 
approximation) is entirely defined by the velocity distribution in the reference 
model: the modelled events therefore cannot adjust to multiple reflections if the 
moveout thereof differs from that of the primary reflections. However, because 
of the stationarities at the zero offset, there is in practice a partial adjustment 
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and the linearized inversion results appear to be contaminated by the presence 
of the correlated noise consisting of the multiple reflections. The seismic data of 
Fig.12A have been modelled by taking account of the multiple reflections (there 
are only 3 primary reflections which reach the zero offset successively at the 
times: 500 ; 1000 ; 1650 ms, the 2 nd arrival having a particularly high amplitude). 

[0086] In the seismic response of the model obtained after conventional 
linearized inversion (see Fig.12B), the high-amplitude multiple arriving at 
approximately 1500 ms at the zero offset has been only weakly attenuated and 
the model is contaminated by the multiple. It is then natural to try to use the 
possibilities afforded by Nemeth's method in order to better discriminate 
multiples and primaries, notably with the improvements provided by our method. 

[0087] We d o scribo h Hereafter the implementation of the method for 1D 
linearized inversion of surface seismic data is described . 

[0088] Experimental protocol 

[0089] It consists in sp e cify i ng specifies : 

the seismic excitation mode: in eu^experiments, the seismic source is a 
source point time-modulated by a function (seismic wavelet) denoted by w(t), 

the various locations where the receivers are positioned: in ewf 
experiments, the pickups are arranged at the depth 10m and cover the offset 
interval [x min =0, x ma x=1500m], with a receiver every Ax=100m ; we-thus obtained 
are seismic traces numbered from 0 to I, 
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the pickups measure the pressure field resulting from the wave 
propagation, 

the recording time:: in ou^-experiments, the pickups measure the 
vibrational state over the time interval [0, T=1 800ms], these data being sampled 
every At=1 ms. W e thus obta i n Thus samples are obtained numbered from 0 to 
N. 

[0090] We call The quantity di n us referred to as the data sample recorded 
for the offset x min + i Ax and at the time n At. W e have of course The relationship : 
T = N At and x ma x = Xmm + I Ax exists . 

[0091] Acquisition of data in accordance with this experimental protocol 

[0092] ©w-eExperiments were carried out from synthetic data: these data 
were obtained by numerical resolution (finite-difference method) of the 2D wave 
equation: 

I d 2 P „<r- 


- V — VP = S(x 7 z-z s )w{t) inmxK * + x[0, T] 


ere dt c 

with the boundary condition: 
P(z=0)=0 

and the initial conditions: 

dP 

P(r = 0) = — (r = 0) = 0 
dt 

an equation wherein: 

x, z, t respectively designate the lateral coordinate, the depth and the 


time, 
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a and c respectively represent the acoustic impedance and propagation 
velocity distributions (functions of the depth). 

[0093] The wavelet selected is a conventional Ricker wavelet of central 
frequency 30 Hz and the impedance a and velocity c distributions are given in 
Figs.13A, 13B. The seismic response of this model therefore cons i sts of js_the 
data represented in Fig.12A. Wo call The quantity F NL is referred to as the 
operator which associates the seismic data with function pair (a(z),c(z)). 

[0094] Specification of the modelling operator 

[0095] W e s e l e ct a A reference model is selected which is described by 
function pair (ct 0 (z),c 0 (z)). In the experiments presented hereafter, w e have 
sel e cted c 0 (z)=c(z) and a 0 (z)=cte (=1) have been selected . W e choos e a s the 
mod e lling op e rator t The Jacobian operator of operator Fnl is chosen as the 
modelling operator at point (<?o(z),Co(z)). In fact, since w e have s ele ct e d 
Co(z)=c(z) has been selected , the only unknown is 8a(z)=a(z)-a 0 (z) and only the 
corresponding component of the Jacobian is important. 

[0096] The modelling operator, which we-eal tis referred to as F(8m), thus is 
the (linear) operator which associates with function pair 6m=(8a(z),6c(z)) the 
samples 5y(Xj,t n ) which are the components of vector F (5m) for Xj=x min + i Ax 
and t n =n At, i=0 ,l ; n=0,...,N. 


31 


612.42948X00 


[0097] Specification of the modelling operator associated with each 
correlated noise 

[0098] We sp e c i fy in th i s paragraph t The procedure is specified herein that 
can be used for modelling each correlated noise. Each correlated noise w ill b e iis 
characterized by correlation directions specific thereto, which we- are assumed 
to be known, at least approximately: these correlation directions are specified 
by means of a field of correlation vectors Cj(x,t) of components (Cj x (x,t), q*(x.t)) 

which has to be defined at any point of the domain [x m in, Xm a x]X[0,T]. W e can, iln 
an equivalent way, sp e c i fy the correlation directions can be specified by means 
of a pencil of correlation lines that wtfJ-represent the field lines associated with 
Cj(x,t). Such correlation lines can be obtained according to the technique 

described in pa t e nts EP - 35 4 ,112 ( US- Patent 4 ,972,383V- filed by the 
applicant assiqnee , from a pick of several phases characteristic of the noise, the 
information being extended to the whole domain [Xmj n , x max ]X[0,T] by means of 
conventional interpolation and/or extrapolation procedures, as described in the 
aforementioned paten t EP - 354,112 f il ed by the applicant . For this second 
application example wherein we — want to — e liminate — elimination of t he 
aforementioned multiple is desired , we p i ck this multiple is picked (the 
amplitude peak for example), which enables us to dof i no defining the arrival time 
variations as a function of the offset and we-defines a correlation line continuum 
by simple vertical translation of the pi ck e d chosen line (under these conditions, 
vector c(x,0 does not depend on t). Specification of the correlation directions 
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amounts to specifying the propagation velocity distribution of the noise Cj(x,t) 


[0099] Unlike the 1 application example described above, wo have here 
noise amplitude variations are along the correlation line. We The assum e 
assumption is that we c a n e st i mate these amplitude variations can be estimated 
from a measurement on the data or by means of theoretical considerations: this 
measurement defines a function g(x). 

[01 00] The noise will thus be modelled by the composition of two operators: 

A transport operator, as in the first application example ; 

An amplitude modulation, jrer that is the multiplication of the solution to 
the transport equation by function g(x). 

[0101] If w e us e again the scheme given for the 1 st application example is used 
(and notably if we use t he same hypotheses and notations are used ), we 
introduc e the following information is introduced : 

a space Bj of noise-generating functions which will be the space of 
support functions pj(t) in [tj min ,tj max ] c [0,T] and that we extend by 0 in the whole 
interval [0,T] ; 

noise modelling operator T,: 

Tj-.fijiOeBj^bjixrfeD 

solution to the transport equation 

Vbj(x 9 t)£j{x,t) = 0 in [x^ 9 x^]x[0,T] 
and meeting the initial conditions: 


that is related to the correlation vector by the relation c.(x,t) = 


c;(m) 



cj(M) 
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bj(x,t = 0) = 0 


[0102] In fact, the transport equation is solved by means of a numerical scheme. 
It is for example possible to use the conventional finite-difference centered 
scheme. W e t Therefore s o l o ct discretization intervals Ax/ and At/ are selected 
(which are not necessarily equal to Ax and At, but which we -are assumed to be 
submultiples of these quantities, even if more general situations can be 
considered) and w e introduce a grid is introduced whose nodes are the points 
of coordinates (x min + i'Ax/, n'At/) with i' and n' e N. The conventional finite- 
difference centered scheme explained hereunder allows, from the initial 
conditions, to calculate step by step the different values taken by function bj(x,t) 
at the various nodes of the grid (to simplify the notations, w e have left out 
subscript j associated with the correlated noise considered has been omitted ). 


[0103] In the above formula, quantity ocr n represents the evaluation of quantity a 
(a is a generic term) at the point of coordinates (Xmm + i'Ax/, n'At/). 


[01 04] Here, discretization intervals Ax/ and At/ have to be selected sufficiently 
small in relation to the wavelength because we want the numerical scheme is 
desired to give a precise approximation of the function solution to the transport 
equation ( w e want no dispersion is desired ). 


2Ax' r n 



■> + b?'+ l - 6J!' +l - if'] + 
+ &S+i - ft?'] =0 
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[01 05] The last stage of the noise modelling procedure consists — m 
s e l e ct i nq selects the samples that belong to the nodes common to the two 
grids and in multiplying them by g(iAx): we thus obtain the various components 
of the vector (of the data space) Tj(Bj) are obtained . 

[0106] Specification of a norm or of a semi-norm in the data space 

[01 07] A first choice would consist in taking t akes as the norm in the data space 
any discrete norm meant to be an approximation (via a quadrature formula) to 
the norm defined on a continuous basis by: 


i 



[01 08] We can also use t The semi-norm can also be used , which is a better 
choice as we shall se e in th e next paragraph discussed below : 



where u is any vector of the data space (i and n are the indices representing the 
sample numbers respectively in space and in time). 

[01 09] Besides, other choices are also possible. 
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[0110]Spec/77caf/on, in each noise-generating functions space, of a norm for 
which each noise modelling operator establishes an isometric relation (or the 
approximation to an isometric relation) between the noise-generating functions 
space and the data space (provided with the semi-norm or defined-m-§§) 

[01 11] We specify i n this paragraph Herein the procedure is specified that can 
be used to define the norm in the noise-generating functions space associated 
with each correlated noise. 

[01 12] A first choice would consist in providing provides each noise-generating 
functions space Bj with any discrete norm meant to be an approximation (via a 
quadrature formula) to the norm defined on a continuous basis by (here again, 
w e have l e ft outt he subscript j has been omitted to simplify the notations): 


OfT.rn.ax \ / *T ^ N 


in which case, using the fact that bj(x,T)=0 for x e [x m jn,x ma x], linear operator Tj 
performs an isometry between Bj and D. However, since calculation of Tjpj is 
carried out numerically, the isometric character of operator Tj will not be perfect 
but only approximate, the approximation being all the better as discretization 
intervals Ax, At, AXj, Atj' are small. The fact that w e th e n hav e only 
approximately an isometry exists will lead to lower performances concerning the 
algorithmic implementation presented in paragraph 8. 
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[01 13] A better choice cons i sts i n prov i d i nq provides space for each noise- 
generating function s spac e function Bj with the semi-norm (here again, w e hav e 
tefroutthe subscript j has been omitted to simplify the notations): 


pii* = 


ax At ^jy(<A*)j \J2 Q (# n+1 + P 1 ? 


[01 14] In this case, operator Tj rigorously performs an isometry between Bj and 
D, at least if the solution to the discretization scheme of the transport equation 
is identically zero for n=N. 

[0115] Definition of the cost function 

[01 16] W e d e fin e t The cost function is defined by the formula as follows: 


C(Sm, ft l9 p 2 ,...fij) = \f nl (m 0 ) + F{Sm) + ]T T.fij - J 
II J* Wd 


[01 17] It can be noted that, in the above function, F NL (m 0 ) is a constant vector. 

[011 8] Seefc/ng the model and the noise-generating functions minimizing the 
cost function, this search being carried out by an algorithmic method taking 
advantage, explicitly or implicitly, of the isometry properties of the noise 
modelling operators (s ee §6) 

[01 19] The algorithmic method uses the specificity of Schur's complement, a 
specificity associated with the isometric character of operators Tj. This can be 
done for example by realizing that minimizing C(6m,pi, P2,..., Pj) amounts to 
minimizing 
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C{Sm) = C(5m,f} x (Sm) 9 P 2 (8m),...fJj(Sm)) 

where (/3 l (Sm),j3 2 (8m),...f3 J (Sm)) minimizes C(5m,p 1( P2,..., Pj) for given 5m. It 

can be noted that the Hessian associated with this quadratic form is Schur's 
complement. C'(5m) can be minimized by means of any optimization method ; 
cost function C'(5m) being quadratic, the use of a conjugate gradient method is 
here particularly appropriate. 

[0120] Determination of (j3 l (m) 9 /3 2 (m),.../3 J (m) is very easy thanks to in view of 
the isometric character of operators Tj, notably if choices (4) and (5) have been 
respectively selected to define and . If there is only one type of noise 

(J=1), which is the case for our illustration, and if the isometry is perfect (which 
is also the case thanks to the choices made), determination of p x {m) simply 
requires evaluation of the gradient of the quadratic form. If there is a 
superposition of several noises (J>1), various algorithms can be considered to 
take advantage of the isometric character of operators Tj as explained for the 
first application example. 

[0121]Figure 14B shows the results obtained (seismic response of the solution 
model) by implementing the method when function g(x) is selected in a very 
simple form (affine function). The improvement in relation to the conventional 
inversion result can be observed (Fig.14A). 

[0122] Implementation examples where the physical parameter whose 
subsurface distribution is modelled is the acoustic impedance have been 
described. It is clear that the method, in its most general definition, can be 
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applied for seeking the distribution of physical quantities affected by correlated 
noises in any heterogeneous medium. 
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CLAIMS 

11) A method of estimating, from data obtained by exploration of a 
zone of a heterogeneous medium, a model representative of a distribution, in 
the zone, of at least one physical quantity, the model being free of a presence 
of correlated noises that may be contained in the data, comprising: 

a) acquiring measurements giving information about physical 
characteristics of the zone by following a predetermined experimental protocol; 

b) specifying a noise modelling operator which associates, with a 
model of each physical quantity, synthetic data that constitute a response of the 
model, the measurements and the synthetic data belonging to a data space; 

c) selecting, for each correlated noise referenced by a subscript j 
ranging from 1 to J, a noise modelling operator which associates a correlated 
noise with a noise-generating function belonging to a predetermined spaced of 
the noise-generating functions (Bj); 

d) specifying a norm or of a semi-norm in the data space; 

e) specifying a semi-norm in the space of the noise-generating 
functions for which each: noise modelling operator establishes substantially an 
isometric relation between the space of noise-generating functions and the data 
space; 

f) defining a cost function quantifying a difference between the 
measurements on one hand and a superposition of a model response and of 
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the correlated noise associated with the noise-generating function on the other 
hand; and 

g) adjustment of the model and of the noise-generating functions by 
minimizing the cost function, by means of an algorithmic method taking 
advantage of isometry of properties of the noise modelling operators. 

12) A method as claimed in claim 1 1 , wherein: 

a distribution as a function of depth of an acoustic impedance in the 
medium is sought, the correlated noises affecting the data are tube waves each 
identified by parameters characterizing their propagation, the measured data 
are VSP data obtained by means of pickups which detect displacement of 
particles in the medium in response to a localized seismic excitation, the 
location of the pickups and a recording time and time sampling points being 
defined, and the selected noise modelling operator associates the synthetic 
data with an acoustic impedance distribution as a function of an evaluated 
depth in traveltime and with vertical stress measured as a function of time at a 
depth of the first pickup. 

13) A method as claimed in claim 12, wherein: 

the cost function quantifying the difference is a square of the semi-norm 
of the difference in the data space. 

14) A method as claimed in claim 12, wherein: 
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adjustment of the model and of the noise-generating functions is 
obtained by means of a block relaxation method for eliminating unknowns 
corresponding to each correlated noise generating function, the block relaxation 
method being implemented within iterations of a quasi-Newtonian algorithm for 
calculation of the model. 

15) A method as claimed in claim 13, wherein: 

adjustment of the model and of the noise-generating functions is 
obtained by means of a block relaxation method for eliminating the unknowns 
corresponding to each correlated noise generating function, this the block 
relaxation method being implemented within the iterations of a quasi-Newtonian 
algorithm for calculation of the model. 

16) A method as claimed in claim 12, wherein: 

numerical calculation of the image of a model by the modelling operator 
is carried out by a numerical solution of a 1D wave equation for the model, by 
selecting values taken by displacement of the particles at locations of pickups 
and at previously defined time sampling points, and by applying an operator 
likely to compensate for spherical divergence and attenuation effects. 

17) A method as claimed in claim 13,wherein: 

numerical calculation of the image of a model by the modelling operator 
is carried out by a numerical solution of the a 1D waves equation for the model 
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considered, by selecting values taken by the displacement of the particles at the 
locations of pickups and at the previously defined time sampling points, and by 
applying an operator likely to compensate for the spherical divergence and 
attenuation effects. 

18) A method as claimed in claim 14, wherein: 

numerical calculation of the image of a model by the modelling operator 
is carried out by a numerical solution of the a 1 D waves equation for the model 
considered, by selecting values taken by the displacement of the particles at the 
locations of pickups and at the previously defined time sampling points, and by 
applying an operator likely to compensate for the spherical divergence and 
attenuation effects. 

19) A method as claimed in claim 15, wherein: 

numerical calculation of the image of a model by the modelling operator 
is carried out by a numerical solution of the a 1D waves equation for the model 
considered, by selecting values taken by the displacement of the particles at the 
locations of pickups and at the previously defined time sampling points, and by 
applying an operator likely to compensate for the spherical divergence and 
attenuation effects. 

20) A method as claimed in claim 12, wherein: 
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the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing a noise transport equation, and a utilized noise 
generating function has an initial condition along an edge of an observation 
zone belonging to the space (Bj) of support time functions in a give time interval. 

21) A method as claimed in claim 13, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 

22) A method as claimed in claim 14, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 

23) A method as claimed in claim 15, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involyed as the function has an initial condition along the an 
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edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval 

24) A method as claimed in claim 16, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 

25) A method as claimed in claim 17, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 

26) A method as claimed in claim 18, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 
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27) A method as claimed in claim 19, wherein: 

the noise modelling operator is a finite-difference centered numerical 
scheme for discretizing the a noise transport equation, and the a utilized noise- 
generating function involved as the function has an initial condition along the an 
edge of the an observation zone belongs belonging to a the space (Bj) 
consisting of the support time functions in a give time interval. 

28) A method as claimed in claim 12, wherein: 
a semi-norm selected for the data space is: 


and a semi-norm selected for the space of the noise-generating function is: 


i 



i 


\m\s = 


Ax At I 



29) A method as claimed in claim 13, wherein: 


the a semi-norm selected for the data space is: 
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and the a semi-norm selected for the space of the noise-generating function 
space is: 



30) A method as claimed in claim 14, wherein: 
the a semi-norm selected for the data space is: 



and the a semi-norm selected for the space of the noise-generating function 
space is: 


i 
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31 ) A method as claimed in claim 1 5, wherein: 
the a semi-norm selected for the data space is: 


Ml* = f Ax A* £ E W +1 + «.")*) 

\ 1=0 n=0 T 2 / 


and the a semi-norm selected for the space of the noise-generating function 
space is: 


H\\ B = (axah ]t + r) 2 ) 


32) A method as claimed in claim 16, wherein: 
the a semi-norm selected for the data space is: 

Hi* = (a* a* J2 E ^tt («r l + o 2 ) 2 

and the a semi-norm selected for the space of the noise-generating function 
space is: 
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33) A method as claimed in claim 17, wherein: 
the a semi-norm selected for the data space is: 

INI* = ( ax a* E E -st « +l + "V 2 ) 2 

and the a semi-norm selected for the space of the noise-generating function 
space is: 


i 



34) A method as claimed in claim 18, wherein: 
the a semi-norm selected for the data space is: 
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and the a semi-norm selected for the space of the noise-generating function 
space is: 


35) A method as claimed in claim 19, wherein: 
the a semi-norm selected for the data space is: 


1Mb = (A* A* £ E « +l + <)*) 

\ i=0 n=0 T 2 / 


and the a semi-norm selected for the space of the noise-generating function 
space is: 
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36) A method as claimed in claim 20, wherein: 
the a semi-norm selected for the data space is: 



and the a semi-norm selected for the space of the noise-generating function 
space is: 



37) A method as claimed in claim 21 , wherein: 
the a semi-norm selected for the data space is: 

1Mb = (Ax At ^fij « +l + 

and the a semi-norm selected for the space of the noise-generating function 
space is: 
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38) A method as claimed in claim 22, wherein: 
the a semi-norm selected for the data space is: 


1Mb = (ax At £ £ -^r « +1 + <) 2 ) 


and the a semi-norm selected for the space of the noise-generating function 
space is: 


PUB - 


39) A method as claimed in claim 23, wherein: 
the a semi-norm selected for the data space is: 
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|M| S = f Ax At £ £ "^T « +l + 

\ i=0 n=0 T 2 / 


and the a semi-norm selected for the space of the noise-generating functions 
space is: 


40) A method as claimed in claim 24, wherein: 
the a semi-norm selected for the data space is: 


1Mb = (ax At £ £ -i^ ( w r +1 + 

\ 2=0 n=0 T 2 / 


and the a semi-norm selected for the space of the noise-generating function 
space is: 
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41 ) A method as claimed in claim 25, wherein: 
the a semi-norm selected for the data space is: 


1Mb = f Ax At £ £ W +1 + <)* 

\ 2=0 n=0 T 2 


and the a semi-norm selected for the space of the noise-generating function 
space is: 


\m\s = ( Ax At i £ +r ) 2 ) 

\ n=0 T 2 / 


42) A method as claimed in claim 26, wherein: 
the a semi-norm selected for the data space is: 



and the a semi-norm selected for the space of the noise-generating function 
space is: 
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43) A method as claimed in claim 27, wherein: 
the a semi-norm selected for the data space is: 


1Mb = ( Ax A* £ E « +1 + 

\ 1=0 n=Q T 2 J 


and the a semi-norm selected for the space of the noise-generating function 
space is: 


P\\B = 


= (a*a*/|;^<*" +1 -m") 2 ) 


44) A method as claimed in claim 1 1 , wherein: 

a distribution of disturbances, in relation to a previously selected 
reference model, of an impedance and of a velocity in the zone of the medium 
is sought, the correlated noises affecting the data are due to multiple reflections 
whose kinematics and amplitude variations with an offset have been previously 
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estimated, the measured data are picked up by seismic surface pickups, the 
location of the pickups, a seismic excitation mode, a recording time and time 
sampling points being defined, and the modelling operator being defined by 
linearization of the waves equation around the model representative of the 
distribution. 

45) A method as claimed in claim 44, wherein: 

the cost function quantifying the difference is the square of the semi- 
norm of the difference in the data space. 

46) A method as claimed in claim 44, wherein: 

adjustment of the model and of the noise-generating functions is 
obtained by a block relaxation method for eliminating unknowns corresponding 
to each correlated noise generating function, the relaxation method being 
implemented within iterations of a conjugate gradient algorithm for calculation of 
the model. 

47) A method as claimed in claim 45, wherein: 

adjustment of the model and of the noise-generating functions is 
obtained by means of a block relaxation method for eliminating the unknowns 
corresponding to each correlated noise generating function, this the relaxation 
method being implemented within the iterations of a conjugate gradient 
algorithm for calculation of the model. 
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ABSTRACT 

[01 23] M e thod A method of forming, from data obtained by exploration of a zone 
of a heterogeneous medium ( a n und e rground zon e for e xampl e ) , a model 
representative of the distribution in the zone of at least one physical quantity 
amply free of the presence of correlated noises that may be contained in the 
data is disclosed . The method essentially comprises the following stages any 
acquisition of data giving information about certain characteristics of the zone, 
specification of a modelling operator which associates the response of the 
model with a model (synthetic data), modelling of each correlated noise by 
applying a specific modelling operator to a noise-generating function (n.g.f.), 
specification of a semi-norm in the data space, specification, in each n.g.f. 
space, of a norm for which each noise modelling operator constitutes an 
isometry, and seeking, by an algorithmic method taking advantage of th o se 
i som e try prop e rti e s, the model and the n.g.f. which minimize a cost function 
quantifying, by means of the semi-norm, the difference between the measured 
data and the superposition of the response of the model and of the noises 
associated with the n.g.f. 

Application: for e xampl e s ee k i ng the distribution, in an und e rground zone, of 
the acoustic imp o d a nc o , of propagat i on v e loc i t ie s, or p e rm e abi l it i es, e tc. 
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